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I. INTRODUCTION 

In the last decade, much interest has been attracted to 
studies of complex networks consisting of dynamical ele- 
ments involved in a set of interactions [H [Ij . Particular 
attention has been paid to problems of synchronization 
in network-organized oscillator systems 0, 3 ■ Synchro- 
nization phenomena are ubiquous in various fields of sci- 
ence and play an important role in the functioning of 
living systems Q ■ Investigations focused on understand- 
ing the relationship between the topological structure of 
a network and its collective synchronous behavior 
Recently, synchronization properties of systems formed 
by phase oscillators on static complex networks, such as 
small- world networks and scale- free networks 0, 01 j 
have been considered. It has also been shown that the 
ability of a network to give rise to synchronous behav- 
ior can be greatly enhanced by exploiting the topologi- 
cal structure emerging from the growth processes [9|, llC| . 
However, full understanding of how the network topology 
affects synchronization of specific dynamical units is still 
an open problem. 

One possible approach is to use evolutionary learn- 
ing mechanisms in order to construct networks with pre- 
scribed dynamical properties. Several models have been 
explored, where dynamical parameters were modified in 
response to the selection pressure via learning algorithms, 
in such a way that the system evolved towards a speci- 
fied goal [nHlg . In our study, this approach is employed 
to design phase oscillator networks with synchronization 
properties. We consider adaptive evolution of a network 
of coupled heterogeneous phase oscillators [13, [13 ■ In 
such systems, heterogeneity of oscillator frequencies com- 
petes with the coupling which favors emergence of coher- 
ent dynamics [1, [13] ■ The question is how to connect 
a set of phase oscillators with given natural frequencies. 
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SO that the resulting network would exhibit the strongest 
synchronization, under the constraint that the total num- 
ber of available links is fixed. 



Previously, a related, but different problem of syn- 
chronization optimization in a network with the fixed 
topology through the modification of connection weigths 
was considered [l9j . Assuming that the system was in 
a phase-locked state, the deterministic steepest descent 
method was used to determine the coupling strengths 
between elements which lead to the best possible phase 
synchronization. In contrast, we consider the systems 
which stay in partially synchronized states (that is, are 
not fully phase-locked) and ask what should be the op- 
timal topology of connections, with each link having the 
same strength. 



To design optimal networks, stochastic Markov Chain 
Monte Carlo (MCMC) method with replica exchange is 
used by us. Large ensembles of optimal networks are 
constructed and their common statistical properties are 
analyzed. As we observe, the typical structure of a 
synchronization-optimized network is strongly dependent 
on its prescribed connectivity. Sparse optimal networks, 
with a small number of links, tend to display a structure 
with relatively high clustering, similar to that found for 
the networks of chaotic maps [20,[1H . As the connectivity 
is increased, synchronization-optimized networks show a 
transition to (approximately) bipartite architectures. 



The paper is organized as follows. In Sec. |TT1 we intro- 
duce a model of heterogeneous phase oscillators occupy- 
ing nodes of a directionally coupled network and define 
the synchonization measure for this system. The opti- 
mization method is also introduced in this section. Con- 
struction of the optimized networks and their statistical 
analysis are performed in Sec. IIIII The results are finally 
discussed in Sec. IIVI 
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II. THE MODEL AND THE OPTIMIZATION 
METHOD 

We consider N oscillators with different natural fre- 
quencies placed onto the nodes of a network. The evolu- 
tion of this system is given by 



dt 



X 



N 



(1) 



where tUi is the natural frequency of oscillator i and A 
is the coupling strength. The weights Uij define the ad- 
jacency matrix a of the interaction network: aij = 1 if 
oscillator i interacts with oscillator j, and a^j = oth- 
erwise. The adjacency matrix is generally asymmetric. 

To quantify synchronization of the oscillators, the Ku- 
ramoto order parameter 



r(t) 



1 

N 



N 



i=l 



(2) 



is employed. Under perfect synchronization, we have 
r = 1, whereas r ^ 0{N~^/'^) in absence of coupling 
for randomly drawn natural frequencies. A second-order 
transition takes place at some critical coupling strength 
A^from the desynchronizcd to the synchronized states 
0. 

To measure the degree of synchronization, we nu- 
merically integrate Eq. ([1]) for given initial conditions 
9i{t = 0) G [0, 27r) and calculate the average modulus 
of r(t) over a long time T, 



i?(a) 




(3) 



where (. . .)init. represents an average over many realiza- 
tions with different initial conditions 0i{O). 

Our aim is to determine the network a which would ex- 
hibit the highest degree of synchronization, provided that 
the total number of links is fixed and a set of natural fre- 
quencies is given. The network construction can be seen 
as an optimization problem. The optimization task is to 
maximize the order parameter and, possibly, bring it to 
unity by changing the network a. An approximate stan- 
dard approach to the problems of complex combinatorial 
optimization, such as the traveling salesman problem, is 
provided by the method of simulated annealing (see, e.g. 

However, we are interested in the statistical prop- 
erties of the synchronization-optimized networks rather 
than in a search for the best-optimized network. If mul- 
tiple samples are generated using conventional optimiza- 
tion methods such as simulated annealing, it is difficult to 
control the probability of the repeated appearance of the 
same (or similar) items in the obtained set of samples. 

To study statistical ensembles of optimized networks, 
the MCMC method p3 - |2^ . which has previously been 
applied to dynamical systems (25143^ . will be used. The 



canonical ensemble average of a network function /(■) is 
introduced as 



/(a)exp(/3j?(a)) 

Zif3) 



(4) 



where Z{f3) = exp(/3i?(a)) is the partition function 
and the parameter (3 plays the role of the inverse tem- 
perature. 

Hence, the problem is reduced to sampling from the 
ensemble with the Gibbs distribution exp(/3i?(a)). Such 
ensemble can be generated, for example, by using the 
Metropolis algorithm [s^, which is the simplest imple- 
mentation of the MCMC method. The Metropolis algo- 
rithm, which we use, is essentially standard. The only 
important difference is that we should simulate the dy- 
namics with a network a at each iterated trial. 

This Metropolis algorithm appears to provide a sim- 
ple and universal way of generating the Gibbs network 
distribution. However, the efficiency of such algorithm 
gets worse when /3 increases, particularly in the case of 
a highly jagged landscape R{a). This deficiency can be 
eliminated by using instead the Replica Exchange Monte 
Carlo (REMC) algorithm, which provides an efficient 
method to investigate systems with rugged free-energy 
landscapes, specifically at low temperatures [33 - l36| . 

In a REMC simulation, a number of replicas {a,„} 
with different inverse temperatures /?„ are evolved in 
parallel. At regular evolution time intervals, the perfor- 
mances of a randomly selected, adjacent pair of repli- 
cas are compared. The running configurations of the 
two selected replicas are exchanged with the probabil- 
ity min [1, exp (A/3Ai?)], where A/3 = Pm+i — Pm is 
the difference of the inverse temperatures of the pair 
and Ai? = R{a.m+i) — i?(wm) is the difference of their 
performances. The exchange of replicas with different 
temperatures effectively imitates repeated heating and 
annealing, thus preventing trapping in the local per- 
formance optima. Note that such stochastic exchange 
algorithm preserves the joint probability distribution 
Hm exp (/3mi?(am))/Z(/3m), so that the unbiased set of 
samples is generated for all inverse temperatures. 

Explicitly, the algorithm is defined as follows: 

1. The states of replicas {a^} are initialized by ran- 
dom networks (which is chosen as a random Erdos 
-Rcnyi network) 

2. The candidate for the next network a^^ at iteration 

(n) 

step n is obtained from the current network a^ by 
rewiring one of its links. A randomly chosen link is 
moved to a randomly chosen link vacancy, so that 
the total number of links remains conserved. 

3. The evolution equations ([1]) for the network a^j are 
integrated using the standard Euler algorithm. The 
order parameter is then calculated and averaged 
over the time interval t G [0, T] and over a fixed 
number of realizations starting from different ran- 
dom initial conditions. Thus, the synchronization 
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FIG. 1: (Color online). The evolution of order parameters of 
coupled oscillator networks during the optimization process. 
The blue solid, red broken, and yellow dotted lines are for 
P = Po, and Pm, respectively. Note that the blue solid line 
is for Po = and, therefore, it corresponds to the networks 
generated by only random rewiring. The parameters are p = 
0.2, iV = 20, A = 1.0, 7 = 0.3, M = 21, 5/3 = 10. 



property i?(aj„) of the candidate network is deter- 
mined. 

4. Next, a random number x G [0,1] is uniformly 
drawn. If 



X < 



exp(/3i?(a^))) 



(n + l) 



the candidate is accepted and taken as a- 
aj^; otherwise nothing is changed, so that a'-"^"'^-' = 

5. At regular evolution time intervals, the perfor- 
mances of a randomly selected, adjacent pair of 
replicas are compared. The running configurations 
of the two selected replicas are exchanged with the 
probability 



1, exp 



{(/3™+i - /3„0(i?(a(:+i)) - i?(w(„"+i))} 



6. Return to Step (2) until the statistical average 
Eq. (HI converges. 



III. NUMERICAL ANALYSIS 

To determine the synchronization degree of a given net- 
work at each iteration step of the optimization procedure, 
equations ([l} were numerically integrated with the time 
increment At = 0.05. Averaging over five independent 
realizations started from different random initial condi- 
tions has been furthermore performed at each iteration 
step. Oscillator ensembles of sizes iV = 10 and 20 were 
considered. Natural frequencies of the oscillators were al- 
ways chosen as Wi — + 2"fi/N , so that they uniformly 
distributed within the interval [—7,7] fioj . 



Initial phases 9i{Q) = 2Tr finit{i) / N uniformly dis- 
tributed inside the interval [0,27r), where fmitii) is a 
random one-to-one mapping between {1, • ■ • , N}. Hence, 
the order parameter at t = always zero. To construct 
initial random networks with a given number K of con- 
nections and, thus, the connectivity p ~ K/N{N — 1), K 
off-diagonal elements of the adjacency matrix were ran- 
domly and independently selected and set equal to unity. 

For time averaging, intervals of length T = 100 and 
200 were typically used. The results did not signifi- 
cantly depend on T when sufficiently large lengths T 
were taken. Using the order parameter, graphs a were 
sampled by the REMC optimization method. In par- 
allel, evolution of M replicas with the inverse temper- 
atures /3m = Sf3 X m, m = 0,1, . . . , M was performed 
(with M = 21 and 6/3 = 10). At each five Monte Carlo 
steps (mcs), the perfomanccs of a randomly chosen pair 
of replicas were compared and exchanged, as described 
above. For display and statistical analysis, sampling at 
each every 50 mcs after a transient of 5000 mcs has been 
undertaken. 



A. Optimization at different temperatures 

Synchronization-optimized networks were obtained by 
running the evolutionary optimization. In this process, 
the order parameter was progressively increasing until a 
stationary state has been achieved. Figure [T] displays 
the optimization processes at different temperatures. As 
clearly seen, when using replicas with the larger inverse 
temperature /3, the larger values of the order parame- 
ter could be reached, although the optimization process 
was then more slow. After the transients, statistical av- 
eraging of the order parameter over the ensemble with 
the Gibbs distribution has been performed, according to 
Eq. Q. _ 

In Fig. [2Ia), the averaged order parameter R is dis- 
played as a function of the connectivity p for several dif- 
ferent inverse temperature f3. The blue solid circle sym- 
bols show the averaged order parameter corresponding 
to the replica with /3o = 0, i.e. for an infinitely high 
temperature. We see that the averaged order parameter 
increases with the network connectivity p even if the net- 
works are produced by only random rewiring. The red 
open circles show the average order parameters for the 
ensemble corresponding to the replicas with the lowest 
inverse temperature /3m- Generally, greater order pa- 
rameters are obtained by running evolution at higher in- 
verse temperatures /3 at any network connectivity p. At 
each connectivity p, the order parameter is gradually in- 
creased with increasing (3 and is approximately saturated 
at f3M- This means that, even if one further increases (3, 
only slight improvements of the averaged order param- 
eter can be expected. Thus, the networks sampled by 
the replica with the largest inverse temperature /3m are 
already yielding a representative optimal ensemble. 

Figure [2I^b) shows, depending on the network conncc- 
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FIG. 2: (Color online). Average order parameters as functions of the network connectivity p. The blue filled circles are for the 
replica j3o, i.e., the ensemble of randomly rewired networks. The red squares, yellow diamonds, green triangles, blue inverted 
triangles, and red open circles are for replicas with /3 = 40, 80, 120, 160, and 200, respectively, (b) Ratio of the average order 
parameters for the synchronization-optimized ensemble with the inverse temperature /3a/ and for /3o = 0. (c) Variance of the 
order parameters. The red squares are for the random rewired ensemble, the blue circles are for the synchronization-optimized 
ensemble. The same parameters as in Fig. [1] 



tivity p, the ratio Rp^j/Rpg of the averaged order param- 
eters sampled by the optimal network ensemble with /3j\/ 
to those obtained for the ensemble with pm'ely random 
rewiring. Since there is no room for the improvement of 
the order parameter when the number of links is small, 
the ratio tends to unity as the connectivity p is decreased. 
On the other hand, when p = 1, global coupling is real- 
ized, for which, under the chosen coupling strength, full 
synchronization occurs. As evidenced by this figure, the 
difference between the synchronization capacities of the 
optimzed and random networks is most pronounced at 
the intermediate connectivities p. 

In Fig. EJc), the mean variance Var[i?]/3 = R'^ — 
of the order parameters at different connectivities p is 
displayed. It can be observed that this mean variance for 
the synchronization-optimized ensemble decreases with 
an increase in the number of links, while the respective 
mean variance for the random rewired ensemble has a 
maximum at p = 0.4. Note that, since the transition 
from the connected to the disconnected random graphs 
occurs at Pc = l/N [1 [|], this behavior is not directly 
related to the topological transition in the network itself. 

To further analyze the behavior of oscillators in 
synchronization-optimized networks, we calculated time- 
averaged frequencies, i.e., winding numbers fli = 
(1/T) 6i(t)dt of all oscillators i. Histograms of dis- 
tributions over the winding numbers were constructed 
by counting the numbers of oscillators with the wind- 
ing number inside a fixed bin interval, Hk = {fli\n6Cl < 
< {k + l)Sn}, where fc = 0, 1, . . . , JsT - 1, A' = 10 
is the number of bins, and 6Cl — 2"f / K is the bin size. 
The winding number as a function of the natural fre- 
quency is shown in Fig. |3ja). The blue circles show 
the entrained cluster with the winding number approxi- 
mately equal to zero. The cluster consists of the elements 
whose natural frequencies are near the mean natural fre- 
quency = 0. While the specific elements of the clus- 
ter and its size depend on a particular network in the 



synchronization-optimized ensemble, there is a statisti- 
cal trend that the entrained cluster consists of the oscil- 
lators in the neighborhood of the zero frequency. This 
is demonstrated by the histogram of winding numbers 
for the synchronization-optimized ensemble in Fig.[3l^b). 
Note that the oscillators arc always ordered according 
to their natural frequencies uji = —7 -I- 2"fi/N which 
monotonously increase with i. We see that all elements 
get divided into two groups, in which w or where 
the winding number is relatively high. For each particu- 
lar network realization, there should be a peak at the fre- 
quency of the entrained cluster. The position of this peak 
depends however on the realization and, as a result, the 
histogram of the winding numbers for the entire ensemble 
shows a broad maximum. This behavior is characteristic 
for relatively low connectivities. The broad peak gradu- 
ally sharpens when the connectivity is increased because 
the size of the cluster increases and fluctuations of the 
winding number become smaller. 



B. Architectures of Synchronization-Optimized 
Networks 



Typical structures of synchronization-optimized net- 
works are shown in Fig. ID When the connectivity p is 
small, such networks usually represent chain fragments. 
At a higher connectivity, the network become more com- 
plexly organized, as shown in Fig.UJ^b). 

To statistically characterize the architecture of con- 
structed networks, ensemble averages of their adjacency 
matrices over the Gibbs ensemble, i.e.. 



a^ = ^aexp(/3i?(a))/Z(/3), 



(5) 



for different connectivities p were computed for /3 ~ , 
as shown in Fig. [S] Clearly, the optimal network struc- 
ture is changing with the number of links. When the 
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FIG. 3: (Color online), (a) Time- averaged frequencies (winding numbers) of the oscillators in a synchronization-optimized 
network, (b) Statistical distribution of winding numbers for the synchronization-optimized ensemble. The parameters are 
p — 0.2, /3 = /3a/; other are the same as in Fig. [T] 
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FIG. 4: (Color online). Two typical realizations of the 
synchronization-optimized network at different connectivities 
(a) p = 0.05 and (b) p = 0.2. The blue (gray) nodes indicate 
entrained oscillators. The numbers in the nodes are indexes of 
the oscillators. The parameters are N — 10, P — /3a/, M — 11, 
and (5/3 = 10. 



number of links is small, the elements of the mean adja- 
cency matrix, obtained by averaging over many realiza- 
tion from the synchronization-optimized ensemble, are 
large near the diagonal. Hence, elements with close nat- 
ural frequency tend to connect and form a chain frag- 
ment. Moreover, oscillators with the natural frequencies 
near the center of the interval are often connected. In- 
creasing the number of links, the network becomes more 
complicated and off-diagonal elements begin to dominate 
instead. The network with the larger p tend to have in- 



terlaced structures, seen in Figs.[5ljb)(c), where the oscil- 
lators with roughly opposite natural frequencies are cou- 
pled. A similar trend towards anti-correlations for the 
oscillators with opposite frequencies has been noticed in 
0, EBl : where a transition from local to global synchro- 
nization under an increase of the coupling strength has 
been obtained using a different optimization method fl6j . 

This structural transition can be understood as fol- 
lows: When connectivity is small, a limited small num- 
ber of available links is better used to connect oscillators 
with frequencies in the middle of the frequency interval, 
where the collective synchronization frequency would lie. 
Indeed, such oscillators can be easily entrained and even 
a single link may be sufficient to synchronize them. If 
connectivity is increased and some further links may be 
used, it would not however be efficient to put them into 
the middle region: the oscillators there are already syn- 
chronized and bringing more connections would not in- 
crease the performance. This means that the additional 
available links should be rather connected to the elements 
in the periphery, outside of the central frequency region. 
If predominantly local connections between the elements 
on each side are established, this would however lead to 
the formation of two clusters, each on a different side 
from the center. Within each cluster, oscillators may get 
synchronized, but oscillations of the two clusters will still 
then be incoherent. Therefore, a better solution would 
consist in establishing pairwisc connections between the 
elements on both sides of the center, i.e. in linking pref- 
erentially the opposite oscillators. This is exactly what 
we observe in Fig. [5] at the higher connectivity p = 0.3. 



C. Degree distributions and cluster organization 

To statistically investigate architectures of designed 
networks, ingoing and outgoing degrees of their nodes 
have been considered and averaged over the ensemble. 
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FIG. 5: (Color online). The upper panels show adjacent matrices averaged over the Gibbs ensemble of synchronization- 
optimized networks [see Eq. ([5])]. The darker color of a matrix element indicates the higher probability of the respective 
connections between the elements. The lower panels display the corresponding network averaged over the Gibbs ensemble. The 
numbers in the circles show indexes of the oscillators. The thickness of the lines connecting the nodes is proportional to the 
frequency of links between them. The network connectivity is (a) p = 0.05, (b) 0.2 and (c) 0.3. Other parameters are same as 
in Fig. [U 



Since the network is colored, i.e, each its node has a dif- 
ferent natural frequency, the mean in- and out-degrec of 
the nodes can be plotted as a function of their natural 
frequency (Fig. [6]). When connectivity p is small, both 
in- and out-degrees averaged over the ensemble have a 
maximum at w = 0, i.e, oscillators having smaller mag- 
nitudes of the natural frequency tend to be mutually 
connected. This unimodal degree distribution is consis- 
tent w^ith the linear chain structure shown in Fig. [5] (a). 
As p is increased, the mean in-degrec distribution be- 
comes bi-modal and oscillators having larger magnitudes 
of the natural frequency tend to have larger out-degrees. 
This tendency becomes stronger when p increases [Fig. [5] 
(b)(c)]. 

Furthermore, we calculated the mean numbers of iso- 
lated nodes as a function of p. The isolated nodes have 
been classificcd into three categories, as those which have 
no in-coming, no out-going, and neither in-coming nor 



out-going connections. The numbers of such isolated 
nodes are, respectively, 

N N 

i=i j=i 

N N 

N N N 

fc=l i=l j = l 

where A{w) is the Kronecker symbol, A{w) = 1 for w = I 
and A{w) = otherwise. We averaged these numbers 
over the Gibbs ensemble for f3 — f3o and (Sm and deter- 
mined the ratio / of the average number of isolated 
nodes in the synchronization-optimized networks to that 
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FIG. 6: (Color online). Degrees, averaged over the the Gibbs ensemble of synchronitaion-optimized networks, as functions of 
the natural frequency of the oscillator. The in-degree fc^^^ , out-degree fc^^^ , and degree k^^^ are plotted by the blue circles, red 
squares and yellow diamonds. A'^ = 10, /3 = /3m; The network connectivity is (a) p = 0.1, (b) p = 0.2 and (c) p = 0.3. Other 
parameters the same as in Fig. [1] 
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FIG. 7: (Color online). The ratio of the number of isolated 
nodes, averaged over the replicas with Pm (synchronization- 
optimized networks) to that averaged over the replicas with /3o 
(randomly rewired networks) as a function of the connectivity 
p. The data for the nodes isolated with respect to incoming 
(blue circles) and outgoing (red squares) connections, as well 
as for the completely isolated nodes (yellow diamonds), is 
shown. The same parameters as in Fig. [T] 



optimized network with larger connectivities may be sim- 
ilar to bipartite graphs. A bipartite graph is a graph 
whose nodes can be divided into two disjoint sets A and 
B, so that every link connects a node in A to a node 
in B and vise versa [37[. To demonstrate that our op- 



timized networks are indeed similar to bipartite graphs, 
we divide all oscillators into two groups A and B with 
the negative and positive natural frequencies. An intra- 
connection is defined as a link between nodes belonging 
to the same group, while an interconnection is a link 
between the nodes in A and B. Thus, the number of 
intraconncctions is given by 



n 



intra 



N/2,N N.N/2 

H-( E + E 

i=lj=JV/2 4=^/2,^=1 

and the number of interconnection is 

N.N 



.winter 



n 



N/2. N/2 

H = ( E + E Kj- 

i=l,i=l i=N /2,j=N /2 



in the networks obtained by random rewiring (see Fig.[7|). 

The results do not depend on the choice of /3m quali- 
tatively. 

When -p is small, the ratio of completely isolated nodes 
is larger than one. This comes from the fact that the links 
are used intensively between the nodes having smaller 
magnitudes of the natural frequency, at the cost of con- 
nections of periphery oscillators. Thus, the number of 
isolated nodes is large. Starting from p ~ 0.23, this ra- 
tio becomes however less than one, so that the optimized 
networks tend to have less completely isolated nodes as 
their random counterparts. We can also notice that the 
relative number of nodes without ingoing connections be- 
comes high at about p ~ 0.23 and then sharply drops 
down. The number of nodes without the outgoing con- 
nections in the optimized networks remains always larger 
than in the random networks. 

As already suggested by Fig.[S](b)(c), synchronization- 



The mean ratio n™*'^''/n™*''° of inter- to intraconnections 
in the synchrony-optimized ensemble for /Ss, /3io, i^is and 
/320 a-s a function of the connectivity p is shown in Fig. [51 
This ratio is smaller than unity when connectivity p is 
small. It increases with p and reaches a maximum in the 
vicinity of the transition point, where the bipartite-like 
structure emerges. Further above the transition point, 
the ratio gradually decreases to unity, since the number of 
links increases until all-to-all connections are established 




D. Closeness, betweenness and clustering 

To characterize network structure quantitatively, we 
calculated the closeness, betweenness and clustering coef- 
ficient Q,!!^- Again, averaging was performed over many 
realizations of synchronization-optimized networks, sam- 
pled with the Gibbs distribution (Eq. 0]). 
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FIG. 8: (Color online). The ratio of inter- to intra- connec- 
tions as a function of the connectivity p for j3 = Ps (blue cir- 
cles), /3io (red squares), /3i5 (yellow diamonds) and /320 (green 
triangles). A'^ — 10. The same other parameters as in Fig. [1] 



The betweenness centrality of a node is the number of 
geodesies (i.e., shortest paths) going through it. If there 
is more than one geodesic between two nodes, the number 
of geodesies which connect these two nodes via a consid- 
ered node is divided by the total number of geodesies 
that connect the two nodes. The betweenness centrality 
is thus defined by 



C" 



s,t s^t 



<Js,t 



where a is the number of shortest paths from node s to 
node t and cr^ t(i;) is the number of shortest paths from 
s to t that pass through node v. 

The closeness centrality of a node specifies how easily 
other nodes can be reached from it (or, in other words, 
how easily it can be reached from the other nodes). It is 
defined as the sum of the lengths of all geodesies leading 
to or from the given node, divided by the total number 
n of nodes minus one. 



where dg{u,v) is geodesic distance between the nodes u 
and V (i.e., the length of the shortest path connecting 
them). 

The clustering coefficient of a node specifies the num- 
ber of neighbours of this node which are in turn mutual 
neighbours. It is defined as 

C*™(i;) = ^, 

where is the degree of a node v and t„ is the number 
of links between its neighbors, Cj" is the number of pairs 
that can be made by using ky neighbors. 

The above properties are defined for each node. To 
characterize the entire network, we average them over all 
nodes. 



In order to quantify differences between 
synchronization-optimized networks and networks gener- 
ated by random rewiring, ratios [w) {w) ^^can 

be used, where is the respective property of net- 
work, such as closeness, betweenness, or clustering, 
is inverse temperature and /3o = 0. In Fig. [HI 
we show these ensemble-averaged network properties 
depending on the connection probability p for several 
inverse temperatures. Obviously, these ratios should 
approach unity at p = or at p = 1, because the 
difference in synchronization of optimized and random 
networks vanishes in these two limits. The ratios for the 
closeness have pronounced minimima in the transition 
region. The ratio in the vicinity of the transition point 
decreases when the performance of optimized network 
increases, i.e., the network ensemble with higher inverse 
temperature. 

On the other hand, the betweenness and clustering co- 
efficient gradually increase with the connectivity p and 
reach a maximum in the transition region. Note that in 
recent work [39j it was found that, both in random and 
scale-free networks, increase the clustering coefficient fa- 
vors formation of oscillator sub-populations synchronized 
at different frequencies. 



IV. CONCLUSIONS 

We have designed synchronization-optimized networks 
with a fixed number of links for a heterogeneous os- 
cillator population. This has been done by using the 
Markov Chain Stochastic Monte Carlo method comple- 
mented by the Replica Exchange algorithm. A transition 
from the linear to bipartite-like networks has been found 
under increasing the number of links. At low connectiv- 
ity, synchronization-optimized networks typically repre- 
sent small chains connecting oscillators with close nat- 
ural frequencies. As the number of links increases, the 
networks become interlaced and oscillators with oppo- 
site natural frequencies tend to be connected. Therefore, 
synchronization-optimized network begin to resemble bi- 
partite graphs. This structural change of synchronion- 
optimized network is clearly revealed through the analy- 
sis of inter- and intraconnections. 

Thus, we have shown that the efficient design of oscil- 
lator networks with the improved synchronization prop- 
erties is possible. The architectures of such optimal net- 
works strongly depend on the constraints, such as the 
total number of links available. Through the appropriate 
rewiring of a network, a strong gain in the synchroniza- 
tion signal can be achieved. 

Although our study has been performed for a sim- 
ple system of phase oscillators, similar evolutionary opti- 
mization methods can be applied to construct networks 
of different origins, where the dynamics of individual os- 
cillators may be significantly more complex. 
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FIG. 9: (Color online). Statistical properties of synchronization-optimized networks as functions of the connectively p. (a) the 
ratio of closeness, averaged over the replicas with j3m to that averaged over the replicas with /3o, (b) the ratio of betweenness, 
and (c) the ratio of transitivity, for different Pm, where m = 10 (blue circles), m = 15 (red squares), and m = 20 (yellow 
diamonds). The same parameters as in Fig. [T] 
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